!***********************************************************
!
program inc_to_tec
!BE CAREFUL WITH FORTRAN COMPILER
!USE OPTION ifort -assume byterecl
!***********************************************************

implicit none



! Variables reading from Incompact3D, real(8) for double precision
integer(4),parameter :: nx=513,ny=193,nz=256,xjicf=372,ibm_y_max=97
real(8), parameter :: lx=24., ly=4., lz=6.

include "declare.f90"
include "tecio.f90"

inquire(file=trim('phim.dat'), exist=iscalar)

if (.not. iscalar) then
  Title="JICF - 3D - no-recycling - no Phi - explicit - averaged"
  Variables="X Y Z u uu uv uw v vv vw w ww"
else
  Title="JICF - 3D - no-recycling - high sc - explicit - averaged"
  Variables="X Y Z t tt u ut uu uv uw v vt vv vw w wt ww"
endif
FName="averaged_JICF.plt"

! Variables for managing files list
nf = ''
nca = command_argument_count()
if (nca==0) then
 write(*,*) 'Programme quits. No files to read / write ...'
 stop
endif

write(*,*) nca,' files to read and export for Tecplot'
allocate (data(nx,ny,nz,nca+3))
allocate (epsilon(nx,ny,nz))
allocate (PassiveVarList(nca+3))
allocate (ValueLocation(nca+3))
allocate (ShareVarFromZone(nca+3))
PassiveVarList=0
ValueLocation=1 !ValueLocation: 0=cell-centered, 1=node-centered
ShareVarFromZone=0
IJKN = nx*ny*nz*(nca+3)

!dx and dz
!dx=lx/nx !periodic
dx=lx/(nx-1) !non-periodic
dz=lz/nz !periodic
!dz=lz/(nz-1) !non-periodic

 
! Reading the mesh
open(21,file='yp.dat',form='formatted')
do j=1,ny
 read(21,*) i,yp(j)
enddo
close(21)
write(*,*) 'Grids y are read'

! reading IBM
if (.false.) then
  OPEN(21,FILE=trim('epsilon.dat'),FORM='UNFORMATTED',&
	ACCESS='DIRECT', RECL=8, STATUS='OLD')!, CONVERT='big_endian')
  data_count = 1
  DO k=1,nz
  DO j=1,ny
  DO i=1,nx
    READ(21,REC=data_count) epsilon(i,j,k)
    data_count = data_count + 1
  ENDDO
  ENDDO
  ENDDO
  CLOSE(21)
else
  DO k=1,nz
  DO j=1,ny
  DO i=1,nx
    tempa=xjicf-i
    tempb=nz/2-k
    r2=dx*dx*real(tempa*tempa)+dz*dz*real(tempb*tempb)
    if ((j .le. ibm_y_max) .and. (r2 > 1.)) then
      epsilon(i,j,k)=1.
    else
      epsilon(i,j,k)=0.
    endif
  ENDDO
  ENDDO
  ENDDO
endif
write(*,*) 'Reading IBM completed'
!epsilon=0.

! Calculating coordinates for Paraview
DO k=1,nz
 DO j=1,ny
  DO i=1,nx
   data(i,j,k,1)=real(   (0+(i-xjicf))*dx) !be careful about (0+(i-1)*step) or (1+(i-1)*step)
   data(i,j,k,2)=real( yp(1+(j-1))   )
   data(i,j,k,3)=real(   (0+(k-nz/2))*dz) !be careful about (0+(k-1)*step) or (1+(k-1)*step)
  ENDDO
 ENDDO
ENDDO

! ATTENTION: SI DATA EN PROVENANCE IDRIS ALORS CONVERT='big_endian' (IF DATA FROM IDRIS THEN CONVERT)

do l=1,nca

 call get_command_argument(l,nf)
 ! reading file
 ! write(*,*) 'Reading file ',trim(nf)
 OPEN(21,FILE=trim(nf),FORM='UNFORMATTED',&
    ACCESS='DIRECT', RECL=8, STATUS='OLD')!, CONVERT='big_endian')
 data_count = 1
 DO k=1,nz
  DO j=1,ny
   DO i=1,nx
    READ(21,REC=data_count) data(i,j,k,l+3)
    if (epsilon(i,j,k)==1.) data(i,j,k,l+3)=0.
    data_count = data_count + 1
   ENDDO
  ENDDO
 ENDDO
 CLOSE(21)
 write(*,*) 'Reading ',trim(nf),' completed'

enddo
! end loop

write(*,*) Variables

err = tecini112 &
(trim(Title)//char(0), &
trim(Variables)//char(0), &
trim(FName)//char(0), &
"."//char(0), &
FileType, &
Debug, &
VIsDouble)

!/*  Ordered Zone */
err = teczne112 &
("Ordered Zone"//char(0), &
ZoneType, &
nx, &
ny, &
nz, &
ICellMax, &
JCellMax, &
KCellMax, &
SolTime, &
StrandID, &
ParentZn, &
IsBlock, &
NFConns, &
FNMode, &
TotalNumFaceNodes, &
TotalNumBndryFaces, &
TotalNumBndryConnections, &
PassiveVarList, &
ValueLocation, &
ShareVarFromZone, &
ShrConn)

err = tecdat112(IJKN, data, DIsDouble)
err = tecend112();

write(*,*) '-----------END-----------'

end program inc_to_tec

